uma introdução
A pesquisa em causalidade busca descobrir uma relação funcional entre duas ou mais variáveis. Por exemplo:
\[ y = f(x) = a + bx \]
No mundo real, observamos essa relação com ruídos. O importante é que esse ruído não tenha relação com outras variáveis que explicam o fenômeno. Exemplo:
\[ y = f(x) = a + bx + \varepsilon, \ \ \varepsilon\sim\mathcal N(0,\sigma^2) \]
Vamos pensar em um problema onde estamos estudando se a grama fica molhada, observando o céu e um regador.
A -> B -> C
. . .
A <- B -> C
. . .
A -> B <- C
Nesse caso, temos de controlar por Idade e sexo
Nesse caso, não é para controlar por pressão
lm() ajusta um modelo linear.anova() serve para testar modelos encaixados.aov() é como se fosse anova() + lm() juntos.tukeyHSD().\[ SQT = SQR + SQE \]
\(SQT\): Soma de quadrados total \(\sum_{ij}(y_{ij}-\bar y)^2\)
\(SQR\): Soma de quadrados dentro de cada grupo \(\sum_{j}\sum_{i}(y_{ij}-\bar y_j)^2\). Quero que isso seja pequeno.
\(SQE\): Soma de quadrados entre grupos \(\sum_{j}(\bar y_{j}-\bar y)^2\). Quero que isso seja grande.
A estatística de teste é dada por
\[ \frac{SQE / gl_E}{SQR / gl_R} \]
aov() e lm()aov() testa o fator todo de uma vez.
lm() testa cada nível do fator individualmente
É possível reproduzir o comportamento de aov() com a utilização de lm(), um modelo nulo e anova()
modelo_lm <- lm(comprimento_bico ~ ilha, data = pinguins)
modelo_lm_nulo <- lm(comprimento_bico ~ 1, data = pinguins)
# Anova serve para comparar modelos
comparacao_modelos <- anova(modelo_lm, modelo_lm_nulo)
comparacao_modelosAnalysis of Variance Table
Model 1: comprimento_bico ~ ilha
Model 2: comprimento_bico ~ 1
Res.Df RSS Df Sum of Sq F Pr(>F)
1 339 8598.6
2 341 10164.2 -2 -1565.6 30.862 4.86e-13 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# A tibble: 2 × 7
term df.residual rss df sumsq statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 comprimento_bico ~ ilha 339 8599. NA NA NA NA
2 comprimento_bico ~ 1 341 10164. -2 -1566. 30.9 4.86e-13
Outros comandos úteis são:
| Função | Descrição |
|---|---|
confint() |
Intervalo de confiança para os parâmetros |
resid() |
Resíduos do modelo |
fitted() |
Valores ajustados |
AIC() |
Critério de informação de Akaike |
model.matrix() |
Matriz de planejamento (matriz X) do modelo |
linearHypotesis() |
Teste de combinações lineares de parâmetros |
vcov() |
Matriz de variância-covariância dos parâmetros |
Objetos de classe formula possuem sintaxe muito conveniente para especificar o modelo estatístico que desejamos ajustar. O símbolo que define esses objetos é o ~.
Estrutura:
Então se o objetivo fosse ajustar o modelo
\[ Y_i = \beta_0 + \beta_1X_i + \varepsilon_i, \]
passaríamos ao R a seguinte fórmula
* e :Utilizamos o símbolo * para introduzir os componentes de interação, além dos componentes aditivos.
Teoricamente teríamos, para Z contínua, o modelo de regressão
\[ Y_i = \beta_0 + \beta_1X_i + \beta_2Z_i + \beta_3X_i*Z_i + \varepsilon_i, \]
ou, para Z categórica, o modelo de ANCOVA
\[ Y_{ij} = \alpha_j + \beta_jX_{ij} + \varepsilon_{ij}, \]
Os operadores aritméticos exercem funções diferentes em fórmulas.
O sinal de + no exemplo induziu em um modelo aditivo em vez de somar X com Z. Para fazer com que eles assumam seus significados aritméticos temos que utilizar a função I().
Exemplo:
Agora sim o componente I(X + Z) representa a soma de X com Z. Outros exemplos: I(X^2), I(log(X + 1)), I(sqrt(X+Z*5)).
| Símbolo | Descrição |
|---|---|
+ X |
inclui a variável X |
- X |
retira a variável X |
X * Z |
inclui X, Z e a interação entre elas |
X : Z |
inclui apenas o componente de interação entre X e Z |
(X + Z + W)^2 |
inclui X, Z, W e as interações 2 a 2 |
I(X + Z) |
Função identidade. Inclui uma variável construída pela soma de X com Z |
X - 1 |
Remove o intercepto (regressão passando pela origem) |
. |
O ponto representa ‘todas as demais variáveis’ |
Frequentemente temos interesse em saber se parâmetros são diferentes de zero ou se são diferentes entre si. Para isto, costumamos efetuar testes do tipo Wald para combinações lineares dos parâmetros.
Para este fim, a função linearHypothesis() do pacote car faz o trabalho.
Outra alternativa é utilizar os testes pareados de Tukey. Nesse caso os valores-p já são ajustados.
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = comprimento_bico ~ ilha, data = pinguins)
$ilha
diff lwr upr p adj
Dream-Biscoe -1.089743 -2.495170 0.3156833 0.1628734
Torgersen-Biscoe -6.306505 -8.203279 -4.4097305 0.0000000
Torgersen-Dream -5.216762 -7.188974 -3.2445487 0.0000000